using System;

using Atomic.Libraries.Mathematics.LinearAlgebra.AlgLib;

namespace Atomic.Libraries.Mathematics.LinearAlgebra
{
	[Serializable]
	public sealed class ComplexEigenvalueDecomposition
	{
		private Vector eigenvalues;
		private ComplexMatrix eigenvectors;

		private ComplexEigenvalueDecomposition()
		{
		}

		/// <summary>
		/// Decomposes a hermitian complex matrix. Only the upper triangular part is used; symmetry is not checked.
		/// </summary>
		public static ComplexEigenvalueDecomposition Decompose(ComplexMatrix a)
		{
			ComplexEigenvalueDecomposition ed = new ComplexEigenvalueDecomposition();
			Decompose(a, out ed.eigenvalues, out ed.eigenvectors);
			return ed;
		}

		/// <summary>
		/// Decomposes a hermitian complex matrix. Only the upper triangular part is used; symmetry is not checked.
		/// </summary>
		public static ComplexEigenvalueDecomposition Decompose(Complex[,] a)
		{
			return Decompose(new ComplexMatrix(a));
		}

		/*public static ComplexEigenvalueDecomposition DecomposeSymmetrized(ComplexMatrix a)
        {
			int n = a.Rows;

			if (a.Columns != n)
			{
				throw new ArgumentException("Matrix isn't a square matrix.");
			}

			//ComplexMatrix
			//return Decompose(0.5 * (a + Matrix.Transpose(a)));
			throw new NotImplementedException();
        }*/

		/// <summary>
		/// Decomposes a hermitian complex matrix. Only the upper triangular part is used; symmetry is not checked. Eigenvalues are sorted in ascending order. Eigenvectors are stored as matrix columns.
		/// </summary>
		public static void Decompose(ComplexMatrix a, out Vector eigenvalues, out ComplexMatrix eigenvectors)
		{
			int n = a.Rows;

			if (a.Columns != n)
			{
				throw new ArgumentException("Matrix isn't a square matrix.");
			}

			AP.Complex[,] a0 = new AP.Complex[n, n];
			for (int i = 0; i < n; i++)
			{
				for (int j = i; j < n; j++)
				{
					a0[i, j] = new AP.Complex(Complex.Re(a[i, j]), Complex.Im(a[i, j]));
				}
			}

			// Eigenvalues are in ascending order.
			double[] d = new double[n];
			AP.Complex[,] z0 = new AP.Complex[n, n];

			evd.hmatrixevd(a0, n, 1, true, ref d, ref z0);

			Complex[,] z = new Complex[n, n];
			for (int i = 0; i < n; i++)
			{
				for (int j = 0; j <= i; j++)
				{
					z[i, j] = new Complex(z0[i, j].x, z0[i, j].y);
				}
			}

			eigenvalues = new Vector(d);
			eigenvectors = new ComplexMatrix(z);
		}

		/// <summary>
		/// Eigenvalues are sorted in ascending order.
		/// </summary>
		public Vector Eigenvalues
		{
			get
			{
				return eigenvalues;
			}
		}

		/// <summary>
		/// Eigenvectors are stored as matrix columns.
		/// </summary>
		public ComplexMatrix Eigenvectors
		{
			get
			{
				return eigenvectors;
			}
		}
	}
}
